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In analyzing systems which depend on uncertain parameters, one tech- 
nique is to partition the uncertain parameter domain into a failure set 
and its complement, and judge the quality of the system by estimating the 
probability of failure. If this is done by a sampling technique such as Monte 
Carlo and the probability of failure is small, accurate approximation can re- 
quire so many sample points that the computational expense is prohibitive. 
Previous work of the authors has shown how to bound the failure event 
by sets of such simple geometry that their probabilities can be calculated 
analytically. In this paper, it is shown how to make use of these failure 
bounding sets and conditional sampling within them to substantially re- 
duce the computational burden of approximating failure probability. It is 
also shown how the use of these sampling techniques improves the confi- 
dence intervals for the failure probability estimate for a given number of 
sample points and how they reduce the number of sample point analyses 
needed to achieve a given level of confidence. 

Keywords: uncertainty propagation, failure probability. 

I. Introduction 

This paper considers a system that depends on parameters subject to uncertainty that 
are modeled probabilistically by modeling the vector of uncertain parameters as a vector 
of random variables. Acceptability of the system is assessed by determining if the system 
satisfies a set of requirements which are prescribed as a set of constraints. 
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The values of the uncertain parameters which result in an unsatisfactory system may 
be agglomerated into an event which will be called the failure domain. A measure of how 
satisfactory the system is can be found in the probability of failure. Computation of the 
exact probability of failure by, for example, computing the integral of the probability density 
function over the failure domain, is usually a computationally intractable problem. This is 
true even if the probability distribution of the uncertain parameter vector is known precisely 
in terms of effectively computable probability density functions or cumulative distribution 
functions, and even if a given realization of the uncertain parameters may be tested for mem- 
bership in the failure domain by an effectively computable procedure (which, in a realistic 
engineering example, often takes a substantial amount of computation). 

A typical technique used to approximate the failure probability is based on sampling the 
uncertain parameter at a large number of points and analyzing the system for each of these 
sample values to determine which sample values fall in the failure domain. The higher the 
precision needed from such an approximation, or the closer the failure probability is to zero, 
the more sample points need to be used; and this technique can also become computationally 
intractable. 

In prior work , 1 the authors have shown how, in this context, to calculate sets whose com- 
plements bound the failure domain and which have such simple geometry (hyper-rectangles 
or hyper-spheres) that their probabilities may be determined analytically for a very useful 
family of possible distributions of the uncertain parameters. These sets have been used to 
find upper bounds for the failure probability. 

In this paper, two techniques are presented which reduce the number of sample points 
which must be subjected to the full computational burden of system analysis while using a 
sampling based method to estimate the failure probability. These are called hybrid methods 
since they both make use of the failure domain bounding sets and also make use of sampling 
techniques. In one, the rejection method , 1 an adequate number of sample points is generated 
in the uncertain parameter space, but a simple test is applied to each to see if it falls outside 
the failure bounding set. If it does, it is unnecessary to go to the computational expense of 
subjecting that point to a full system analysis; i.e., that point is rejected for analysis. In 
the other, the conditional sampling method, the only sample points generated fall within the 
failure bounding set. These points are distributed according to the portion of the distribution 
of the uncertain parameters which is supported by the failure bounding set. This reduces the 
number of points needed (substantially, if the event bounding the failure domain itself has 
small probability) to achieve a given level of confidence in the failure probability estimate 
because, in contrast to the rejection based method, no sample points are ever generated 
that need to be rejected. The conditional distribution technique can be considered to be 
an importance sampling technique 2-4 where the distribution used for sampling is determined 
using information obtained from a failure bounding event (set) which has been found by 
solving an optimization problem, and from which samples can be drawn analytically. 

It is then shown how using these sampling techniques can reduce the number of samples 
needed (vs. full space sampling) to achieve a given level of confidence or, conversely, how 
more confidence can be achieved for an estimate based on a given number of sample points. 
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II. Basic Notions 


The objective of this paper is the analysis of systems having a parametric mathematical 
model. These parameters, denoted by the vector p, are considered to be uncertain. This 
uncertainty is modeled by considering the vector of uncertain parameters to be a vector 
of random variables whose distribution is specified by a joint probability density function 
(PDF) f p (p) whose support set, also called the uncertainty set, is denoted by A p . Hereafter, 
the terms “uncertainty set” and “support set” will be used interchangeably. By specifying 
the uncertainty model, it is implied that A p contains the true value(s) of p. In principle, the 
uncertainty set can be unbounded and can even be the entire parameter space. In practical 
applications, the choice of this model is usually made by a discipline expert. However, 
the theory presented herein acknowledges that such a choice may be fairly arbitrary, and 
addresses the need for quantifying the level of tolerance of the system to uncertainty in 
p. Any member of the uncertainty set is called a realization. The nominal parameter value, 
denoted by p, is a parameter realization regarded as a good deterministic representation of p. 
Requirements on the system are prescribed by means of the collection of constraint functions, 
g(p) < 0 (this and all vector inequalities are considered satisfied if they are satisfied in every 
component). 

A system is considered robust if the constraints are satisfied for all p € A p . For systems 
which are non- robust, the failure probability to be approximated subsequently is a metric of 
the degree of non- robustness. 

The Failure Domain T is given by a 


[jr j , 

3 = 1 

where the scalar components of the vector g(p) are (g i(p) , g 2 (p) , • • • ,g q (p)), and 

F 3 = {P : 9j(p) > 0}- 

Here, is the failure domain corresponding to the jth requirement, and T is the failure 
domain for all requirements. Usually, the equation g - = 0 partitions the parameter space into 
two parts, the failure domain corresponding to the jth requirement, and its complement 
set, the Non- failure Domain corresponding to the jth requirement, (fF 3 ) c (the superscript c 
denotes the set complement operator). 

III. Set Deformations 

Failure set bounding is accomplished by expanding or contracting a specially chosen 
Reference Set. The Reference Set will be denoted by H and will be chosen as a hyper-sphere 
or hyper-rectangle centered at the nominal parameter p. The hyper-sphere of radius R 
centered at p, denoted by S(p, R ), is defined by 

S(p,R) = { p ■ \\p-p\\ < R}- 

a Throughout this paper, super-indices are used to denote a particular vector or set while sub-indices refer 
to vector components; e.g., p\ is the zth component of the vector p-b 
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The hyper- rectangle centered at p whose shape is determined by the vector m > 0 of the 
half-lengths of the sides is denoted R(p, m) and is defined by 

R(p, m) = {p : Pl e [pi - m u p i + mj, 1 < i < s} 

where s = dirn(p). Choice of these reference sets is motivated by the analyst’s perception 
of relative uncertainty of the parameters. Hyper-spherical sets consider parameters with 
similar levels of uncertainty. They could also be used for parameters with dissimilar levels 
of uncertainty if some scaling is used. However, whether scaling is used or not, a degree of 
dependence is introduced. For instance, the range of variability of one parameter depends 
on the values taken on by the others. Rectangular sets permit the consideration of dissimilar 
levels of uncertainty without the need for scaling and introduce no such dependence. Such 
levels are attained by making m, proportional to the level of uncertainty in p t . 

In previous work, 1 a measure of robustness has been assigned to the system based on 
measuring how much the reference set can be deformed before intersecting the failure domain. 
In that study, deformations of are limited to Homothetic Deformations, by which is meant 
that is expanded or contracted preserving proportionality and orientation, while p remains 
fixed. Thus, if the hyper-sphere S( p, R) is deformed by a factor of a, the resulting homothet 
is S(p, aR ), while if the hyper-rectangle R{p, m) is deformed by a factor of a, the resulting 
homothet is R(p,am). Expansions are accomplished when a > 1, while contractions result 
when a < 1. 

Throughout this paper, two uncertainty sets will be called Proportional if one of the two 
sets can be formed from the other by an homothetic deformation. Call the factor of this 
deformation the Similitude Ratio and denote it by a. 

In what follows, it will be assumed that the nominal parameter value is feasible; i.e., 
g (p) < 0 . Intuitively, one imagines that a set proportional to the reference set is being de- 
formed with respect to the nominal parameter point until its boundary touches the boundary 
of the failure domain; i.e., until at least one member of the deformed set is on the verge of 
violating one or more of the requirements in g. Any point where the deforming set touches 
the failure domain is a Critical Parameter Value (CPV). A CPV, which will be denoted 
hereafter by p, might not be unique. The expanded set is called the Maximal Feasible Set 
(MFS) and it is denoted by A4 P . The Critical Similitude Ratio , denoted by a, is the Simili- 
tude Ratio of that deformation, and the parametric safety margin, to be defined later, is a 
dimensional metric that quantifies the size of the MFS. 

Formulations that enable the deformation of hyper-spherical and hyper-rectangular ref- 
erence sets are introduced next. 

III. A. Deformation of hyper-spheres 


The CPV for the jth requirement is given by 


p 1 = argmin { \\p - p\\ : g^p) > 0} , 

p 

(1) 

while the overall CPV is given by 


p = p\ 

(2) 

where 


i = argmin [\\fd — p\ } . 

(3) 
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Hence, the CPV problem is solved for each individual requirement, and, from these, the 
overall CPV is selected which is closest to the nominal parameter point. Additional con- 
straints can be imposed on Equation (1) to prevent the CPV from assuming values leading 
to infeasible systems; e.g., if p x is a mass, the constraint p 1 > 0 must be added. 

The CPV is the p value on the boundary of the failure domain closest to the nominal 
parameter point. The Spherical Parametric Safety Margin (PSM), denoted by ps, is ||p— p||, 
and the Spherical Maximal Feasible Set (MFS) is Ad = S(p,ps). Ad is the largest sphere 
centered at p which does not intersect the failure domain. 

III.B. Deformation of hyper-rectangles 

Recall that the infinity norm in a finite dimensional space is defined as H^H 00 = sup^ { | a? j | } . 
Let us define the m-scaled infinity norm as ||a;||“ = supdlahl/uz.*}. A distance between the 
vectors x and y can be defined as || a? — y ||“. Using this distance, the “unit ball” centered 
at p is just R(p, m). 

The CPV corresponding to the deformation of R(p,m) results from using Equations (1- 
3) after replacing the Euclidean norm with the m-scaled infinity norm defined above. Hence, 
the CPV is solved for each individual constraint function, and the overall CPV is selected 
to be the individual CPV which is closest to the nominal parameter point according to the 
m-scaled infinity norm. For instance, using the definition of the m-scaled infinity norm in 
Equation (1) leads to 


f nr\ 

p 1 = argmin < max — — : qAp) > 0 

p \i <k<s m k “ 

The “max” can be eliminated from the objective function, thereby eliminating potential 
derivative discontinuities, by introducing the Similitude Ratio a defined earlier: 

(p j , a?) = argmin{a : gAp ) > 0, \p t — p { \ < am h 1 < i < s} (4) 

p,a 

As before, additional constraints can be imposed on Equation (4) to prevent the CPV from 
assuming values leading to infeasible systems. The interpretation of p given in the previous 
section can be naturally extended to the rectangular case. Figure 1 shows a sketch with 
relevant metrics. 

The Rectangular Parametric Safety Margin (PSM), denoted by p R , is p R = a||m||, and 
the Rectangular Maximal Feasible Set (MFS) is A A — R (p, am). M. is the largest homothet 
of R(p, m ) which does not intersect the failure domain. 

IV. Failure Probability: Upper Bounds and Hybrid Methods for 

Estimation 

In general, the estimation of failure probabilities is difficult since it requires evaluating 
a multi-dimensional integral over a complex integration domain. Since the direct compu- 
tation of P[F] is usually impractical, one needs to use approximations such as First Order 
Reliability Method (FORM) 2 or Monte Carlo sampling. 5 Sampling techniques tend to be 
computationally intensive. 
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Recall that a probabilistic uncertainty model is prescribed by the joint PDF f p (p), whose 
support set is A p ; or, equivalently, by the joint cumulative distribution function (CDF) 
F p (p). Strategies for the calculation of upper bounds to the probability of failure based on 
the deformations of sets in p-space and -u-space (to be defined) are presented next. 

IV. A. Upper bounds to Failure Probability 

Let 

ijj = P[M C ] = 1 - P[M\, (5) 

where M. is the maximal feasible set resulting from a homothetic deformation of a reference 
set fh This number 0 is an upper bound to P[T}. 

When the components of p are independent random variables and Q is hyper-rectangular, 
0 can be calculated by a simple formula 1 which requires using evaluations at only two points 
of each of the CDFs of the individual components of p. 

One important special case of parametric uncertainty arises when the parameter vector 
has the standard normal distribution. This is not only a useful model in itself, but is 
also used in analyzing other parameter uncertainty models when a probability preserving 
transformation U can be found so that, when u = U(p), then f u is the standard normal 
distribution. Under such a transformation, if T is the failure set in p-space, U{T) is the 
failure set in the corresponding u-space, and the constraint function g in p-space transforms 
into the constraint function G(u) = p(U _1 (u)). Such a transformation is used in many 
standard techniques such as FORM. In this paper, the term u-space refers to an uncertain 
parameter space in which the uncertainty is the standard normal distribution. Note that if 
a system has both p-space and w-space representations connected by a transformation U, 
reference set choices are made independently in the two spaces and are not connected through 
the transformation U, since U will usually fail to preserve the geometric form needed for 
the failure probability bounding calculations made in this paper. This also means different 
bounds may be achieved in the two representations, and no general rule is known to determine 
in advance which is better. 
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IV. B. Hybrid Methods for Estimating Failure Probability 

The developments that follow are applicable to both p - and u-spaces. 

In Monte Carlo sampling, a representative sample {p 1 ,p 2 , ■ ■ ■ ,p N } of the random vector 
p is generated, the constraint function values g(p l ), i — 1, . . . N, are calculated, and P[T] 
is estimated as the fraction of these constraint vector values which place p l in the failure 
domain. The value of N is dictated by the confidence desired in the resulting approximation 
(the higher the confidence, the bigger N must be) and P[F] (the closer P{J~\ is to 0, the 
bigger N must be). Since P[F] is presumably unknown in advance, it might be necessary to 
pick a trial value for N, generate N samples, and estimate P[F] and calculate the confidence 
interval at some level of confidence (frequently, 95%) for this estimate, and decide if more 
samples must be analyzed to gain an acceptable level of confidence. (As an aside, note that 
the same problem arises if P[P] is close to 1. However, as a practical matter, it is not 
anticipated that failure probability estimates will be made in cases where failure is almost 
certain.) 

Very small probabilities arise in practical engineering examples. For example, the FAA 
Federal Aviation Regulations (FARS, 14 CFR) Part 25 - Airworthiness Standards: Trans- 
port category airplanes, Section 1309 - Equipment, systems, and installations includes the 
specification that 


(b) The airplane systems and associated components, considered separately and 
in relation to other systems, must be designed so that - 

(1) The occurrence of any failure condition which would prevent the continued 
safe flight and landing of the airplane is extremely improbable . . . 

In Advisory Circular No. 25. 1309-1 A, the FAA quantifies their phrase “Extremely Improb- 
able” as “having a probability on the order of 10 -9 or less.” 

Hybrid Methods, to be introduced subsequently, approximate the failure probability by 
combining sampling methods with the failure set bounds. Applications of the rejection based 
method and a conditional sampling method to problems in w-space have been made previ- 
ously. 1 This paper extends these ideas to problems in p-space, and presents improved and 
extended versions of the conditional sampling techniques. The present p-space conditional 
sampling techniques can be applied in rt-space. Using the present techniques in rt-space 
constitutes an improvement over the methods given in the Ref. 1. 

One advantage of doing the analysis in p-space is that, as long as the constraint function 
g and the reference set 0 have not changed, the p-space MFS, denoted here by Xi p , remains 
the same even if the uncertainty distribution f p assigned to p changes. However, if the u- 
space has been defined as the image of p-space under a probability preserving transformation 
U and a rt-space MFS, denoted here by M. u , has been calculated using the derived constraint 
function G(u) = g{U~ l {u )), changing f p changes U and consequently G, so a new MFS 
M. u must be calculated. 

Sampling based methods depend on the generation of sample sets from specified dis- 
tributions. One general technique is reviewed here. Suppose that X = (Ad, X 2 , . . . , X s ) 
is a vector of s mutually independent random variables whose CDF is F(x \,x< 2 , ■ ■ ■ ,x s ) = 
F 1 (xi)F 2 (x 2 ) • • • F s (x s ). It is desired to generate a sample set of N vectors from the distribu- 
tion F. First, generate a sample {u 1 , v 2 , . . . , v N } distributed uniformly over an s-dimensional 

7 of 17 


American Institute of Aeronautics and Astronautics 



hypercube which is the Cartesian product of s copies of the interval (0, 1). There is an abun- 
dance of software to find such samples either using pseudo-random number generators (as 
would be done in Monte Carlo analysis) or using deterministic schemes such as Hammers- 
ley Sequence Sampling or Halton Sequence Leaped 6 for generating so-called low discrepancy 
sequences. Then, if x{ = F~ l {vl) } the set { x 1 ,x 2 , . . . ,x N } is such a sample set. Standard 
libraries for statistical analysis, such as the Matlab® Statistics Toolbox, contain software to 
evaluate the PDF, CDF, and inverse CDF of many commonly used distributions. 

IV.B.l. Rejection Method 

The simplest of the hybrid methods is the rejection method. In the rejection method, a 
p-spaee sample {p l ,p 2 , ■ ■ ■ ,p N } is first generated according to the distribution f p (p). 

Once the MFS A4 has been determined, an arbitrary sample p‘ can be easily tested for 
membership in this MFS. If p l e Af, it is unnecessary to go to the computational expense 
of evaluating g at p\ since the outcome is now known. For hyper-rectangular reference sets, 
where JA = R(p,am), the membership test for p l E M, is that every component of the 
vector | p l —p\ be bounded by the corresponding component of dm; i.e., || p l — p||“ < a. The 
only situation in which we use hyper-spherical reference sets is in mspace with its standard 
normal distribution. Consideration is limited to hyper-spheres centered at the origin. In this 
case, M. = S u (0,aR), and the membership test for u l E M, is ||w*|| < aR. 

For ease of discussion, it is assumed that the sample points in {p l jP 2 , ■ ■ ■ ,p N } have been 
renumbered so that {p l , . . . ,p n } C M c and {p n+l , . . . ,p N } C M. Since T C M c , and 
recalling the definition of in Equation (5), the conditional probability equation P[A \ B } = 
P[A n B]/P[B] implies (with A = T and B = Ai c ) that 

P[P\ =j)P[T | M c \. 

P[P | M c ] is estimated by ^2'- =1 R (p 1 ) /n where X is the indicator function of T which is 
given by: 

t ( tj ) = / 1 if P e A 
\ 0 otherwise 

Notice that determining which value R(p) takes on requires evaluating g(p). This provides 
the estimate: 

/ n 

( 6 ) 

' l -i 
1=1 

This sum involves only the members of the sample set which are in A4 C . Thus, it can be 
considered that the points of the original sample which fall in A4 have been rejected. This 
motivates the use of the name Rejection Method. Notice also that, since can be calculated 
directly from knowledge of A4 and f p , and since n/N is an estimate of jj, the expected value 
of n, the number of sample points falling inside M. c , can be controlled by the choice of N. 

Equation (6) provides an approximation to P[P\ when is calculated using the overall 
CPV p, and to P[T when is calculated using the individual CPV p J . Note that P[P\ > 
maxlPfX 7 ]}. The membership test above allows identifying the n samples out of N to be 
evaluated by the indicator function. Since the MFS is used in this fashion, this variant of the 
hybrid method is called rejection-based. Other variants 7 of the hybrid method are available. 

Since this method avoids evaluating g at all samples within the MFS, its efficiency in- 
creases with the magnitude of the parametric safety margin. The accuracy of the hybrid 
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method using n samples is comparable to the accuracy of the Monte Carlo method using 
N = [n/'ijj] samples, where [x] is the integer for which [x] < x < [x] + 1. Hence, the larger 
the MFS, the larger the numerical gain. 

IV. B. 2. Conditional Sampling 

The rejection method can save a lot of computational effort on system analysis (evaluation 
of the constraint functions). However, if the failure probability is very small, it is still 
necessary to generate a very large number of samples (the N of the previous section) in 
order to have a sufficiently large number (the n of the previous section) of samples falling 
in T to provide a good estimate of P[P\. The same considerations for determining N for 
Monte Carlo sampling apply to determining N for the rejection method. This would require 
an unmanageable number of samples if P[tF\ were on the order of, say, ICC 9 as in the FAA 
specification mentioned earlier. 

The next technique presented is Conditional Sampling. In this technique, only sample 
points outside the MFS Ad are generated. These are generated according to the portion of the 
distribution of the uncertain parameters which has Ad as its support (see Equation (7)). In 
this way, the expense of generating N — n points which end up falling within Ad and of testing 
them for membership in Ad is avoided. The computational cost of generating each sample 
point using the proposed conditional sampling techniques for either the hyper-rectangular or 
hyper-spherical Ad can be expected to be larger than the cost of generating a sample point 
unconditionally from the distribution of the uncertain parameter. So conditional sampling 
becomes preferable to the rejection method when P[Ad] is sufficiently close to 1. 

Versions are given for MFSs which are hyper-rectangles and MFSs which are hyper- 
spheres. The version for hyper-rectangles requires only that the components of the random 
vector be independent random variables, and is stated for p-space. It could be applied to 
■u-space since the random variables there also have independent components. The hyper- 
spherical version applies only to w-space hyper-spheres centered at the origin. 

If components of a random vector p are independent random variables, then samples 
of the vector can be generated by sampling each component independently. The problem 
to be overcome with the conditional sampling techniques given here is that the conditional 
distributions used no longer describe independent random variables. 

The two conditional sampling techniques given here are examples of what is called in 
the literature importance sampling. See Refs. [2-4] and their bibliographies. In importance 
sampling, the uncertain parameters are sampled, not according to the distribution / which 
actually describes their uncertainty, but according to a different distribution h chosen (typi- 
cally) to give more samples in areas important to the analysis. In order to make the resulting 
analysis valid for the distribution /, the contribution of each sample point x generated ac- 
cording to the h distribution is weighted by the correction factor /(x)/h(x). For this to be a 
valid procedure, care must be exercised in the choice of h ; for example, consideration must 
be given to the support set of / and to the probabilistic estimate being accomplished by the 
sampling procedure. In the present case, the h distributions are the conditional distributions 
used for picking the conditional samples, and the ratio /(x)/h(x) is constant over the range 
of x values of interest, and that constant is the reciprocal of the probability of some bounding 
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set for the failure domain. Specifically, if _Ad c is a failure bounding set and if} = P[Ad c ]: 


h(x) 


0 x E M 

x^M 


( 7 ) 


Some authors have referred to this special case of importance sampling as “Partition of 
the Region.” 


Hyper-rectangles: 

Recall that s has been defined to be dim(p). Here it is assumed that the components 
(p 1 , . . . ,p s ) of the random vector p of uncertain parameters are independent random vari- 
ables. This means that the joint PDF of p and the joint CDF of p are products of the 
respective marginal functions: 


fpip) = f pl {pl) • 

" f Ps {Ps ) 

(8) 

F p(p) = F Pl (Pi) ■ 

■■FpsiPs) 

(9) 


The MFS is a hyper-rectangle: 


M = R(p, am) 


For purposes of this section, it is more convenient to represent A"! as a Cartesian product of 
intervals: 

M. — Ii x • • • x I s 


where 

Ij [o. j , bj \ , j 1, • • • , s 

with a,j = pj — amj and b 3 = p j + arrij . 

The conditional PDF to be sampled from is 


h(p) = f P {p | M c ), 


(we write h instead of h p to reduce notational clutter subsequently when subscripts are added 
to represent various marginal distributions) and is represented by the equation: 


I 0 pe M 

h{P) = 4#t otherwise 

l P[M C } 

Denote the corresponding CDF by H(-). 

A sample point from this conditional distribution is generated one coordinate at a time, 
with coordinate s being generated first, then s — 1, and so on until a value has been generated 
for each coordinate. To do this, it is necessary to refer to the “partial hyper-rectangles” 


M k = h x • • • x I k 


so that M. s = M.. 


10 of 17 


American Institute of Aeronautics and Astronautics 



Recursive formulas give the marginal probabilities of these partial hyper-rectangles and 
their complements. This means that probabilities of M. k and Ai k are given with respect to 
the marginal distribution in fc- dimensional space of the partial random vector (p ly . . . ,p k ). 
The basis step of these recursive formulas is: 


P[M{\ = F Pl (b l )-F Pl (a l ) 

P[ :M;] = fp,(a I ) + (l-F Pl (6,)) 

For k > 1, the recursive step is: 

P\M k ] = (F Pt (b 0 - F n (a k )) 

P[Ml\ = P\M%_ i] + P[M k - J (F Ph (a k ) + (1 - F n (b„))) 

This form of the formulas is chosen rather than calculating P[M. c k ] as 1 — P[M. k \ to avoid 
the numerical pitfalls of subtracting two nearly equal numbers which could arise if the MFS 
had probability nearly equal to 1. Of course, the same problem can arise in calculations of 
the form 1 — F{b). Any relief from numerical problems of this form must depend on specific 
properties of the individual CDFs. For example, if F is the standard normal distribution, 
1 — F(b) can be calculated as F(—b) to avoid the problem. Using software for the standard 
normal CDF from the Matlab® Statistics Toolbox on a 64 bit computer with ANSI standard 
floating point arithmetic for F and b = 8.25, F(—b ) is computed to at least 12 digits accuracy 
while the computation represented by 1 — F (b) is in error by about 40%. 

The sample construction process also uses the marginal distributions of h on each of the 
spaces formed by the first k components of p. So, if 

fkipi, ■ • • , pj = f Pl Gpi) ■■■f Pa ( Pk), 

this marginal conditional PDF is 

h k (Pi, ■ ■ ■ ,P k ) = fk(Pi, ■ ■ ■ ,P k I M c k ), 
and is represented by the equations: 

( 0 p G M k 

WPl ’-- P ‘ )= otherwise 

l P\Ml\ 

Denote the corresponding CDF by H k (p l , . . . ,p k ). 

The generation of a given coordinate of a sample vector depends on where the previously 
generated coordinates lie. Suppose that the coordinates (p k+l , ■ ■ ■ ,p s ) have been generated. 
If any one of the coordinates Pj $ Ij, k + 1 < j < s, no matter how the remaining coordinates 
are picked, the sample point will fall in A4 C , as desired. This means that sample value p k 
may be picked by sampling from the distribution F p . If, however, all of the coordinates 
Pj G Ij, k + 1 < j < s (and this case applies vacuously in case k = s), the sample must be 
chosen according to a distribution which apportions the proper weight to whether the sample 
point p k falls within or outside I k . This requires sampling from H kyk which is the marginal 
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distribution of H k for its last coordinate, h ktk (p k ) is calculated by integrating h k (p 1 , ■ ■ ■ ,p k ) 
over the half-space 

{{Pi, ■■■,Pk)\ ~ oo < Pi < oo, ■ ■ ■ j -oo < p k _ k < OO, -oo <p k < p* k }. 


This produces the formula: 


H k ,k(Pk ) 


-^Pfc(Pfc) 

P[M c k ] 

p Pfc (K) - nMj 


K- < ak 


a k < Pi < h 


h < Pi 


By using P Pfc \ H kk may be evaluated as needed, and the required sample generated as noted 
previously. 

As previously remarked, this technique can be applied in rt-space for arbitrary hyper- 
rectangles. 


Origin centered hyper-spheres in n-space: 


The problem considered here is to sample from the s-dirriensional standard normal dis- 
tribution conditional on the sample points falling outside an origin centered hyper-sphere in 
■u-spaee having radius R. One sampling technique for this case has been presented in earlier 
work. 1 The present work is an improvement over the previous in that sample points are not 
confined to fall on a limited number of lines through the origin, so are distributed in a more 
representative fashion. The key to this is making use of the CDF A s of the random variable 
||ii||. This is known as the ^-distribution with s degrees of freedom. A s is given 1 (see also 
Ref. [8], p. 168, for the corresponding ^-distribution CDF formula) by: 


A s(r) = { 


erf (72) - \Jl ((Syn + (Syi! + • • • + fn) e 7-2/2 if r> 0 ,s odd 
1 “ ((Syn + (SyiT + ' ' ' + fi + x ) e_r2/2 if r > °, s even 

0 otherwise 


Here, n\\ is the double factorial, defined as: 


f n ' 

(n - 2) ■ 

• • 5 • 3 • 1 

n > 0 and odd 

( n ■ 

1 n - 2) • 

• • 6 • 4 • 2 

n > 0 and even 

l 1 



n = —1, 0 


The conditional distribution of the magnitude of the s-dimensional standard normal 
random vector u subject to the condition that r = ||u|| > R is given by: 


A sA r ) 


0 

A 8 (r)-A 3 (R) 
1-A S (R) 


r < R 
r > R 
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To generate n samples of the standard normal random vector of dimension s, first gen- 
erate n unit vectors uniformly distributed over the surface of the unit sphere. One way to 
accomplish that is to first generate n vectors {iq, . . . , v n } in s-dimensional Euclidean space 
which are sampled from the standard normal distribution. The rows of the matrix generated 
by the Matlab® command v = randn(n,s) ; can be used here. Then normalize them by 
setting Wi = Vj/||t;j||. Then sample n numbers {ri, . . . ,r n } from the distribution A The 
desired conditional sample is then calculated as iq = ViWi , 1 < i < n. 

To calculate {ri,...,r n }, one can start with numbers {aq, . . . , sampled from the 
uniform distribution on the real interval (0, 1), and solve each of the (non-linear) equations 
A s,R.( r i) — for r i- This can be accomplished using available software (for example, in the 

Matlab® Statistics Toolbox) for the inverse of the y 2 CDF. 

IV.B.3. Confidence Intervals 

A few remarks on the numerical features of the hybrid method are now in order. If P is a 
sampling-based estimate of P[P], the corresponding 100(1 — a)% confidence interval [L, U ], 
satisfies 

P[L < P < U] — 1 — a. 

Confidence intervals relate P with P[P\ by bounding the true, unknown value of the ap- 
proximation with defined certainty. This technique is illustrated in the literature. 9 Here, 
confidence interval studies for traditional Monte Carlo are compared with those for the hybrid 
method. A model problem is chosen with a two dimensional vector of uncertain parameters 
having a standard normal distribution; i.e., a two dimensional w-space. There is a single 
constraint function which is affine; i.e., linear plus a constant. With this simple model, the 
probability of failure can be determined analytically, and can be adjusted to any desired 
level by adjusting the coefficients in the constraint function. Since P is a binomial random 
variable, we can solve for L and U after applying the binomial test. 

For the first study, the constraint function coefficients are chosen so that the MFS M. 
is an origin centered circle (2- dimensional hypersphere) of radius is 3.25. This dictates that 
P[M\ = 0.9949, P[M C ] = 5.086 x 10“ 3 , and P[P] = 5.72 x 10 -4 . Figure 2 illustrates how the 
width of the confidence interval varies with the required number of g evaluations, namely n, 
for this problem, using both Monte Carlo sampling and the hybrid method. 

The computational improvements of the hybrid method are apparent. The improvement 
in efficiency, denoted by 5n, is given by the horizontal distance between the two curves for a 
fixed U — L value. 5n at a given U — L value is the additional number of samples a Monte 
Carlo simulation requires to get the same confidence interval width. The improvement in 
accuracy, denoted by 5a, is given by the vertical distance between the curves for a fixed 
number of g evaluations. 5a at a given n value is how much wider the confidence interval of 
the Monte Carlo estimate is as compared to the confidence interval of the hybrid method. 
For the particular case shown, improvements of about two orders of magnitude in 5n and 5a 
are attained. 

Metrics for evaluating the relative improvements in efficiency and accuracy, denoted by 
Tj n and Tj a hereafter, are introduced next. For a fixed U — L value, the number of Monte 
Carlo samples at which the constraint function must be evaluated is r) n times as large as the 
corresponding number for the hybrid method. On the other hand, for a given value of n, 
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Figure 2. Number of g evaluations vs. width of the 95% confidence interval for P[!F] = 5.72 x 10 4 . 


the width of the confidence interval of the Monte Carlo estimate is rj a times as large as the 
width of the confidence interval for the Hybrid Method. 

Figure 3 illustrates the dependence of r/ n and rj a on P[P}. In the calculation of rj n , we 
use a = 0.05 and U — L = P[P]/ 10 for P[P\ < 0.5, and U — L = (1 — P[P])/10 otherwise. 
Hence, we set the width of 95% confidence interval to be 10% of the failure probability. In 
the calculation of g a we assume n = 1000. It can be seen that the numerical gain of the 
hybrid method for both rj n and rj a increases with the distance of P[P\ from 0.5. While g a 
reaches 6 when P[T] = 0.005, r] n grows even faster reaching 16.8 at P[P\ = 0.01. In general, 
the numerical improvements are proportional to the size of the MFS. Such improvements 
are substantial when P[P] is relatively large or relatively small, the latter being the targeted 
range in practical, safety-driven applications. When A4 = 0, which is the worst case, the 
hybrid method is equivalent to Monte Carlo. In the last two figures we have purposely left 
out the number of g evaluations required to find the CPV since this number depends on 
the initial condition of the numerical search. For this representative problem the number 
of g evaluations to find the CPV ranged from 9 to 31 and averaged 15. As compared to 
the number of function evaluations due to sampling, this contribution to n is practically 
negligible. Even though the same trends have been observed in large dimensional problems, 
a formal analysis of the method is beyond the scope of this paper. 

V. Example 

In order to be able to evaluate the results of a conditional sampling numerical experi- 
ment, a simple model is chosen which allows analytic determination of P[P] ■ The uncertain 
parameter vector is a 4-dimensional standard normal random vector. A single constraint 
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?[F] 


Figure 3. P[P] vs. r] n (dashed line) and ija (solid line). 


function is used, G(u) = (a, u) — b where (a, u) denotes the vector inner product of a and 
u , a is chosen to be (10,1,1,1), and b = f3\\a\\ where /3 = 5.9978. These choices make 
P[P] = 10~ 9 . This probability is so small that it is not practical to estimate using full space 
sampling methods such as the hybrid rejection method or the Monte Carlo method. For 
example, it would take 1,442,695,000 random samples to have a 50% chance of finding even 
one failure point. To estimate this P{J~] with any significance would take tens or hundreds 
of billions of random samples. This is not currently practical. 

However, P[T\ can be estimated using the spherical MFS version of the conditional 
sampling technique with reasonable precision with much fewer sample points. In this case, 
M. = S u (0,/3), so P[Ai c ] = 293.0 x 10~ 9 . Thus, it would be expected that about one in 
every 293 conditional samples is a failure point. The failure probability can be determined 
with a manageable number of sample points. 

The conditional sampling estimate of P[P] was made using a conditional sample of 
100,000 points. A total of 100 independent repetitions of the sampling estimate were made. 
A plot, shown in Fig. 4, of the relative errors between the estimates and the true failure 
probability shows that the Conditional Sampling technique has captured the failure proba- 
bility to a reasonable level. For example, 92 of the 100 trials show a relative error within the 
±10% range and the worst relative error is about 14%. 

Confidence interval analyses at the 95% confidence level for each of these 100 trials shows 
that the true failure probability falls within the confidence interval in 93 of the trials. 

On a desktop system running Matlab® 7.4.0.336 (R2007a) on an Intel® Pentium® 4 CPU 
3.40GHz computer, it took about 6 minutes to generate and analyze the 10 million sample 
points. One must note that the analysis in this simple example is computationally negligible. 
It is also of interest to note that, if all 10 million samples are analyzed as a single sample set, 
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Figure 4. Relative error in Conditional Sampling estimates of P[JF]. 


the resulting estimate of failure probability is 0.9926 x 10 -9 for a relative error of -0.74%. 
The 95% confidence interval for the 10 million point sample is [0.982 x ICC 9 , 1.003 x ICC 9 ]. 
Note that, to get equivalent results to this 10 million point conditional sampling analysis 
using the Monte Carlo method or even the hybrid rejection method would require about 
34.13 trillion sample points. Based on generating and analyzing 10 million points every 6 
minutes, this would take about 39 years. 

VI. Summary 

Techniques from the literature have been reviewed for using optimization methods to find 
upper bounds of failure probabilities by bounding the failure set (constraint violation set) by 
a geometrically simple set whose probability can be calculated analytically. We have shown 
how to use that failure bounding set to reduce the computational burden for estimating 
failure probability by sampling the uncertain parameters. Some of these techniques may be 
used in p-space while previous methods were limited to w-space. Improvements have been 
made to previous w-space methods. The rejection method generates a full uncertainty set 
sample, but rejects any points which fall outside the failure domain bounding set before 
subjecting the remaining points to the analysis necessary to determine of they are constraint 
violation points. Conditional sampling techniques are given in which only sample points 
within the failure domain bounding set are generated, saving the expense of generating 
points outside the failure domain bounding set. Probability estimates based on these samples 
constitute an application of importance sampling. It is shown how the use of either technique 
based on failure domain bounding improves the ratio of confidence in the failure probability 
approximation per sample point subject to full system analysis. Finally, an example shows 
that Conditional Sampling can be enabling technology in case of need to estimate extremely 
small probabilities. 
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